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Abstract We propose a finite volume method on general meshes for the discretization of 
a degenerate parabolic convection-reaction-diffusion equation. Equations of this type arise in 
many contexts, such as the modeling of contaminant transport in porous media. We discretize 
the diffusion term, which can be anisotropic and heterogeneous, via a hybrid finite volume 
scheme. We construct a partially upwind scheme for the convection term. We consider a wide 
range of unstructured possibly non-matching polygonal meshes in arbitrary space dimension. 
The only assumption on the mesh is that the volume elements must be star-shaped. The 
scheme is fully implicit in time, it is locally conservative and robust with respect to the Peclet 
number. We obtain a convergence result based upon a priori estimates and the Frechet- 
Kolmogorov compactness theorem. 



In this paper we study a finite volume method on general meshes for degenerate parabolic 
convection-reaction-diffusion equations of the form 



Equations of this type arise in particular in the modeling of contaminant transport in ground- 
water. The unknown function u represents the concentration of the species, which diffuses 
and is transported by the groundwater. An essential element in our study is the processus 
of adsorbtion by a porous skeleton, which is supposed to be very fast. More particularly we 
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1 Introduction 




(i) 



suppose that the dissolved and the absorbed parts of the species are in equilibrium; this is 
modeled by the function (3, where 0' may be infinite in several points. The matrix A is a 
possibly anisotropic and heterogeneous diffusion-dispersion tensor, V is the velocity field, the 
function F stands for the chemical reactions, and q is the source term. We suppose that the 
mesh is quite general, and possibly nonmatching. Therefore, also in view of the anisotropy in 
the diffusion term, we can not apply the standard finite volume method [T3] , 

Finite volume schemes have often been applied to the equation (CQ), see e.g. [2], [5], [T5]. The 
upwind discretization of the convection term permits finite volume schemes to be stable in 
convection dominated case, however standard finite volume schemes do not permit to handle 
anisotropic diffusion on general meshes. On the other hand finite element method allows a 
very simple discretization of full diffusion tensors, they were used a lot for the discretization of 
equation ([1]), see e.g. [5J, [9], [TU]- A possible solution is to split equation ([[]) into a hyperbolic 
part and a parabolic part, by means of an operator splitting method; one can find such an 
analysis in [21], [22], where the advection term was treated by the method of characteristics. 
The other quite intuitive idea is to take "best from both worlds" [17J, which leads to combined 
finite volume-finite element schemes; we refer to [17] for this approach. In order to solve this 
class of equations, Eymard, Hilhorst and Vohralfk |T7] discretize the diffusion term by means 
of piecewise linear nonconforming (Crouzeix-Raviart) finite elements over a triangularization 
of the space domain, or using the stiffness matrix of the hybridization of the lowest order 
Raviart-Thomas mixed finite element method. The other terms are discretized by means 
of a finite volume scheme on a dual mesh, where the dual volumes are constructed around 
the sides of the original triangularization. In the second paper of Eymard et al. [18J the time 
evolution, convection, reaction, and sources terms are discretized on a given grid, which can be 
nonmatching and can contain nonconvex elements, by means of a cell-centered finite volume 
method. In order to discretize the diffusion term, they construct a conforming simplicial mesh 
with vertices given by the original grid and use the finite element method. In this way, the 
scheme is fully consistent and the discrete solution is naturally continuous across the interfaces 
between the subdomains with nonmatching grids, without introducing supplementary equations 
and unknowns nor interpolating the discrete solutions at the interfaces. 

The finite volume methods for the discretization of anisotropic diffusion on general meshes is a 
subject of wide interest (see for instance the results of the benchmarks organized at the FVCA 
5 conference [FVCA5]). We refer to [T2J for a detailed analysis of three recently developed 
families of schemes, namely the Mimetic Finite Difference scheme, the Hybrid Finite Volume 
scheme and the Mixed Finite Volume, which turn out to be quite similar. The most important 
feature of these methods is their accurate approximation of anisotropic diffusion even in highly 
heterogenous cases. In this paper, we apply a recent method based upon the finite volume 
method on general meshes developed by Eymard, Gallouet et Herbin [13J, whereas we use 
a slightly modified upwind scheme for the approximation of the convection term. The time 
discretization is based upon a completely implicit finite difference scheme. 

The organization of this paper is as follows. We describe the numerical scheme in Section 2. 
We show the existence and uniqueness of the solution of the discrete scheme and prove a priori 
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estimates for the discrete solution in L°°(0, T; L 2 (Q)) and in a discrete space analogous to the 
space L 2 (0, T; in Section 3. In Section 4, we prove an estimate on differences of time 

translates whereas we establish an estimate on differences of space translates in Section 5. 
These estimates imply a relative compactness property of sequences of approximate solutions by 
the Frechet-Kolmogorov theorem. We deduce the strong convergence in L 2 of the approximate 
solutions to the unique solution of the continuous problem in Section 6. For the proofs, we 
apply methods inspired upon those of [13] and [14] . In Section 7, we finally present results 
of numerical tests, which confirm the validity of the numerical method. 



2 The numerical scheme 

We consider the parabolic degenerate convection-diffusion-reaction problem 

( d/3(u) 



(P){ 



()j V • (A(x) Vu) + V • (V(x)«) + F(u) = g(x, t), (x, t) G Q T , 
u(x,t) =0 x G dQ, te (0,T), 



u(x, 0) = tto(x), x G Q, 

where Q is a bounded open connected polyhedral subset of W 1 , d G N \ {0}, T > and 
Q T = f2 x (0, T). We suppose that the following hypotheses are satisfied: 

(Hi) (3 G C(R), (3(0) = is a strictly increasing function, which satisfies the growth condition 
\(3(a) - f3(b)\ ^ 0\a — b\, (3 > for all a,b el; moreover there exist P > and Cp, such 
that \/3(u)\ ^ Cp for |m| ^ P and /3 is a Lipschitz continuous with a constant (3 for \u\ ^ P; 
("H 2 ) A is a measurable function from Vt to JAdQSL), where A / l ( i(R) denotes the set of d X d 
symmetric matrices, such that for a.e. x G Vt the set of its eigenvalues is included in [A m , Am], 
where A m , Am G are such that < A ^ A m (x) ^ Am(x) ^ A; 

(n 3 ) V G #(div, fi) fl L°°(fi) is such that V • V ^ a.e. in Q; 
(Hi) Uq G 

(^ 5 ) F G C(R), F(0) = and there exists M > such that uF(u) > and F(u) is 
Lipschitz continuous with constant L F for all u < or u > M; moreover we suppose that F 
does not decrease too fast i.e. there exists F > such that (F(u)—F(v))(u—v) ^ —F_(u—v) 2 
for all w,u6l; 
(He) q G L 2 (Q T ). 

We now present a definition of a weak solution of Problem (V). 

Definition 2.1. We say that a function u is a weak solution of Problem (V) if 

(i) u G L 2 (0,T; Hq (£l)); 

(ii) f3(u)eL°°(0,T;L 2 ({l)); 

(Hi) u satisfies the integral equality 

/3(u)(f t dxdt- [ /3(u )(p(-, 0) dx+ [ [ AVu-Vip dxdt 
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— / / uV ■ Vy2 dx.dt + / F(u)<p dxdt — I q<p dx.dt 
Jo Jn Jo Jn Jo Jn 

for all if G L 2 (0,T;H%(n)) with <p t G L°°(Q T ) } <p(-,T) = 0. 

Remark 2.1. In the case that the reaction function F is nondecreasing, the uniqueness of 
the weak solution of Problem (V) follows from [24]. 

In order to describe the numerical scheme we introduce below some notations related to the 
space and time discretization. 

Definition 2.2. (Space discretization) Let VL be a polyhedral open bounded connected 
subset ofR d , with d G N\ {0}, and dVl = Q\Q its boundary. A discretization ofVt, denoted 
by V, is defined as the triplet V = (At, 8, V), where: 

1. At is a finite family of non empty connex open disjoint subsets ofVt (the "control volumes") 
such that Vt = [J K€M K. For any K G M, let dK = K\K be the boundary of K; we define 
m(K) > as the measure of K and Uk ss the diameter of K. 

2. £ is a finite family of disjoint subsets of Vt (the "edges" of the mesh), such that, for 
all a G 8, a is a non empty open subset of a hyperplane ofM. d , whose (d — l)-dimensional 
measure m(a) is strictly positive. We also assume that, for all K G At, there exists a subset 
8k of 8 such that dK = {J a££K o. For each o G 8, we set Ai a = {K G Ai\o G 8k}- 
We then assume that, for all a £ 8, either Ai a has exactly one element and then a G <9f2 
(the set of these interfaces called boundary interfaces, is denoted by 8 ext ) or JW a has exactly 
two elements (the set of these interfaces called interior interfaces, is denoted by 8 int ). For all 
a G 8, we denote by x CT the barycenter of a. For all K G JM. and a G 8k, we denote by n^ i(T 
the outward normal unit vector. 

3. V is a family of points of Vt indexed by Ai, denoted by V = (~&k)k&m< sucn that for all 
K G Ai, Xk £ K; moreover K is assumed to be x K -star-shaped, which means that for all 
x G K, there holds [x K , x] G K. Denoting by d K ^ the Euclidean distance between x K and 
the hyperplane containing a, one assumes that dK,a > 0. We denote by Dk,v the cone of 
vertex x^ and basis a. 

Next we introduce some extra notations related to the mesh. The size of the discretization V 
is defined by 

h v = sup diam(K); (2) 
KeM 

moreover we define 

a s dK,a h K s. 

&x> = max max , max ). (3) 

ae£ in t,{K,L}=M < r d L:(7 KeMa,tre£ K d K ,a 

Thus imposing a uniform bound on 9d forces the meshes to be sufficiently regular. As it was 
done in [T4] we associate with the mesh the following spaces of discrete unknowns 



X v = {((v k )k&m, (va)aes),v K eR,v a G R} 
Xv,o = {v G X v such that (v a ) ae£ext = 0}. 



(4) 
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Moreover, for each function (p = <p(x) smooth enough we define P-pip G Xp in following way 

(Pv<p)k = <p(*k) for all K G M, 
(Pvv)a = V?(x CT ) for all a e £. 

Definition 2.3. (Time discretization) We divide the time interval (0, T) into N equal time 
steps of length St = T/N such that 

6t<0/F, (5) 
where St is the uniform time step defined by St — t n — £ n -i- 

Remark 2.2. For the sake of simplicity, we restrict our study to the case of constant time 
steps. Nevertheless all results presented below can be easily extended to the case of a non 
uniform time discretization. 

After formally integrating the first equation of (V) on the domain K x (i n _ 1; i n ) for each 
K G Ai and n — 1, . . . , N, we obtain 

f /3(w(x, t n )) - P(u(x, f 1 " 1 )) dx + f f (-AVw + Vu) ■ n Kt(T d^dt 

rt n r rt n r 

+ / F(u) dxdt = q dxdt. 

For all K G Ai and all a G £k we define Vk<j — I V ■ n K ^dj and q't = — - / [ q dyidt 

Ja ' St m{K) J tn _ x J K 

We use an upwind scheme in order to approximate the convective term, since it can possi- 
bly dominate the diffusion term; the diffusive flux — / AVw • wk^I is approximated by 

J a 

a function of the form F ^ a {u n ^ , where u n = {{u%)KeMA u V)^s), ar| d where the numeri- 
cal flux FK, a (u n ) is defined by formula ( 125]) below. The time implicit finite volume scheme 
corresponding to Problem (V) is given by: 

(i) The initial condition 



u k = — } T s\ / M o(x) dx, (6) 
Ik 



m{K) 



for all K eM. 

(ii) The discrete equations 



m{K){/3{u n K ) - /3K" 1 )) + St +St^ V ^ u k* 

+St m(K)F(u n K ) = St m(K)q n K , 



(7) 



for all K G Ai. Unlike in the case of the standard upwind scheme, we define the value u\ a 
as follows. For all K G Ai and a G £k we set 

u n K , if V K ,a > 



Zn — —J Ki " r K-,cr ^ u (o\ 

' CT ~ u" if PV, < o. w 
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We also define 

V+„=~(V Kia + \V K , a \) and V^ = \{V K , a -\V K><T \) (9) 

which lead to 

V K ^ = Vl a u n K + V^K (10) 

The definition of flBD seems natural since we also take the unknowns associated with the mesh 
faces. It has an important advantage that the unknowns in the equation ([7]) are associated 
with a single control volume (see Remark I7~I|) : moreover the numerical experiments presented 
in Section 7 show that the upwind scheme (|5]) also preserves the approximate solution from 
unphysical oscillations in the convection dominated case. Finally, we remark that for each 
time step the number of equations is card(A4), whereas the number of discrete unknowns is 
equal to card(Ai) + card(£). Therefore we need to introduce card(E) additional equations 
corresponding to the interface values. For boundary faces these equations are obtained by 
writing the discrete analog of the Dirichlet boundary condition 

(iii) < = for all aeS ext . (11) 

For interior faces, we follow the main idea of the finite volume method by imposing the local 
conservation of the discrete fluxes 



(iv) (F Ka (u n ) + V K>a u^ a ) + (F L>a (u n ) + V L!<T ul >a ) = (12) 

for all a £ £ int with M. a = {K, L}. We will define below Fk, in some more detail, but we 
first give an alternative variational formulation of the discrete scheme (i)-(iv). Let {t> n } ne N 
be an arbitrary sequence of elements of X© ; multiplying equation ([7]) by v K and summing 
on all control volumes KeM leads to: 

KeM ' ' KeM ae£ K 

+ m(K)v n K F(u n K ) = ™(K)v n K q n K . 

KeM KeM 

Using ( Tl2~]) , we obtain that 

<( F l<Au n ) + V K ,atfkJ = f0r aN V U e X V, (13) 

KeM oe£ K 

which yields the following discrete weak formulation: 
Let u Q K be defined by: 

u° K = ——— [ Wx) rfx for all K eM (14) 

m{K) J K 
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For each n G {1, . . . , N} find u n G Xx>,o sue/? that for all v n G Xd ]0 : 

KeM dt (15) 



+ £ m(K)v n K F(u n K ) = m(K)v n K q n K , 



KeM KeM 



with 



and 



<v,u> F = ^2 ^2 ( vk ~ v °) f kA u ) 

KeM oe£ K 



<v,u > T = ^ ^2( v k~ v a )V K>a u KjCr . (17) 
KeM ae£ K 

Remark that the problem (i)-(iv) is equivalent to (TT5i) - (TT71) . Indeed, let 5ij be the Kroneker 
symbol, by setting — for all a G £, and v' K = Skk 1 for all K' G Ai and for a given K 
one recover (ii), and setting v K = for all K G Ai and v a i = 5 aa i for all a' E £ yields (iv). 
The homogeneous Dirichlet boundary condition (iii) follows from the fact that u 11 G X v o . In 
order to complete the numerical scheme we still have to express the discrete flux Fk,c m terms 
of the discrete unknowns. For this purpose we use the SUSHI scheme proposed in [14|: the 
idea is based upon the identification of the numerical fluxes Fk,o through the mesh dependent 
bilinear form, using the expression of a discrete gradient. We first define 

V*ru = — ^2 - VAT G Ai, VmG X v . ( 18 ) 

m ^ > ae£ K 

Remark that the geometrical relation 

^ m(o-)n A> (x (T - yi K ) T = m(K)Id (19) 

ae£x 

holds for each K. Let (p(x) be a function, piecewise linear on the control volumes of the mesh. 
In view of ( TI9~]) one has V kPv{^p) = Vy^x)^.^- We also remark that 

m{o)n K , a = ^2 nA > ^ 7 = / V1 rfx = °> 

which means that the coefficient of uk in ffl8j) is equal to zero; thus, a reconstruction of the 
discrete gradient solely based on ffIB]) cannot lead to a coercive discrete bilinear form in the 
general case. Therefore we introduce the additional term 

V K ,aU = V K u + R K ,aU ■ n X)tT , (20) 

where 

Rk, u = - — {u<r - u K - V K u ■ (x ff -x x )), (21) 

U>K.<r 
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for some oik > 0, which should be chosen in a suitable way. If we choose ax = Vd for all 
K E Ai in the simple case that A is a scalar and that the mesh that satisfies the orthogonality 

property n Kcr = — — — — , we obtain the usual two point scheme. Nevertheless, it may be 

useful to optimize the choice of ax as it is done in |3J. We then define the discrete gradient 
Vv u as the piecewise constant function equal to V^ jCr u in the cone D K(T with vertex x K and 
basis a 

Note that the term Rr,u is a second order error term, which vanishes for piecewise linear 
functions. Moreover, the relation (Tl^j) together with (l2~Tj) implies that 

rn(D K>(T )R Kj<7 (u)n Kt(T = for all K E M and for all u E X v , (22) 

cre£ K 

which in turn implies that 

/ Vx>u c?x = m{K)V kU. 
Jk 

The discrete gradient defined above satisfies the following strong consistency property. 



Lemma 2.1. Let V be a discretization ofVt in sense of Definition 12.21 moreover let 6 ^ 8x> 
be given. Then for all (p E C 2 (VL), there exist a positive constant C only depending on d, 8 
and ip such that 

\\V v P V (f - V^H^oc^d ^ Ch v . 

The proof of this Lemma is given in [13]. The numerical flux is implicitly defined by the 
relation 

< v,u > F = ^2 ( Vk ~ v v) f kA u ) = / Vp?; • A(x)V©m dx. (23) 
It can also be defined explicitly; in order to do so, we write the discrete gradient in the form 

Vk,„u= -MA')y CT<T ', (24) 

where y aa ' is defined by 

m(o-) Vd m(a) , 

——w K . a + - — (1 r—rn K ^ ■ (x CT - XK))n K ,a if a = a , 

m{K) ' dK,a m{K) 

Vd m[u) 

n K ,*< ~ i TTpr^-Ky ■ (x CT - x A ')n^ jCT otherwise. 

m(K) ' dK,arn{K) 

We obtain that for all u,v E Xx> 

I V©u(x) • A(x)V^(x) = S A< k( u «- u k)( v °'- v k), 

•' n KeM <re£ K a'e£ K 
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with 

A k = ^'^ ■ A A>"y CT ' V and A A> „ = I A(x) rfx. 

<r"ee K k ,o n 

The local matrices (A™ ) a a'ee K are symmetric, and the numerical flux is then defined by 

F K ,a( u ) = Yl A k'^k -u a ,). (25) 
a'e£ K 

Next we prove some useful properties of the mesh depending bilinear forms introduced pre- 
viously. In particular we prove that < •, • > F is continuous and coercive, and that < •, • > T 
is continuous and nonnegative. The space X v defined in fl4]) is equipped with the following 
semi-norm. 

Definition 2.4. Let V = (Ai,£,V) be a discretization of Vt in the sense of Definition \2.2\ 
then for all v E Xx> we define 

\v\x= E E^K-^ (26) 



which is a norm on the space X V q. Let us also define a discrete analog of || ■ norm. 

Definition 2.5. (The discrete space H M {yt)) Let 1 ^ p < oo and let V = (M,S,V) 
be a discretization of Vt in the sense of Definition \2.2[ Let H M (Q) C L 2 (fi) be the set of 
piecewise constant functions on the control volumes of the mesh M. for each v e Hm (Vl) we 
define v K = w(x)| xeA -. 

For all v e Hm{&) 3nd for all a 6 £i n t with M. a = {K, L} we define D a v = \vk — Vl\ snd 
d a = dx,a + di,a, 3r >d for all o G S ext with M. a = {K}, we set D a v = \vk\ snd d a = dx,a- 
We then define the following family of norms 

\\<p,M =EE m(a)d K>a (^) p ; (27) 

KeM <reS K ° 



so that in particular 



\1,2,M - E m ( a ) d KA'-J- 



Next we recall two results from [14j which we will use below. The following lemma shows the 
equivalence between the semi-norm in Xx> and the L 2 -norm of the discrete gradient. 

Lemma 2.2. Let V be a discretization ofVt in the sense of Definition ^. 21 and let 9 ^ 9-p be 
given. Then there exists C\ > and C 2 > only depending on 9 and d such that 

Ci\v\ x < || V-dV 1 1 L2(n) ^ C 2 \v\ x for all v e X v . 
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Lemma 2.3. LetV be a discretization ofVt in the sense of Definition ^. 21 then there holds 

\\v\\l,2,M < Mx for a// w e -^X>,0- 

Next we show that the bilinear forms defined in (TIB]) and (TT7|) satisfy continuity and coercivity 
properties. 

Lemma 2.4. /.efD 6e a discretization ofVt in the sense of Definition ^. 21 anc/ /et 6 1 ^ #£> 6e 
given, then: 

(i) There exist positive constants C\ and a which do not depend on h such that 

I < U,V >p I ^ CillilxM 



\x 



and 



< U,U >f^ a \ u \x 



for all u, v G X-u- 

(ii) There exist a positive constant C 2 which does not depend on h that 

I < U, V > T I ^ C 2 |«|xkU 

and 

< U,U >t^ 

for all u,u£ X V fl. 

Proof, (i) Using the definition of the numerical flux (I2"3"j) and in view of ("H 2 ) and Lemma I2~2l 

\<u,v> F \ = \ I V v u ■ \(x)V v v\ dx ^ X\\V v u\\ L 2 {n) \\V v v\\ L 2 {n) < Ci\u\ x \v\ x ; 
Jn 

on the other hand we have that 

<u,u> F = / Vj)U ■ A(x)Vx>u ^ X d^\\V-r>u\\ 2 L 2^ ^ C 2 \u\ x . 
Jn 

(ii) By the definition (TT7|) we have that 

< u> v >T = ^ ^2( v k- v a )V Kt<T W£-. 

Using the definition ([8]) one can write: 

< u ,v > T = ^2 ^2 V kA v k ~ v a )u K + ^2 ^2 V kA v k - v a )u a , 

KeMa££ K ,V K „^0 KeM creS K ,VV CT <0 
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which implies 

< U,V > T = V K) a(vK - V a )u K ^ ^ V^ j(T (f^ - f CT )(M X - (28) 

Using the Cauchy-Schwarz inequality and the bound d^ )(T ^ hx> we have that 

I < «, " >t Sv^.||V|U» (I1) (£ £ mW <!«p£)i(£ ^!!!M^. 4) i 

and 

| < u,v > T | < ||V|| i o (n)(Vd- Mx|M|l2(c) + diam(Vt) ■ \v\ x \u\ x ) 
since cf^m^) = m{K)d. In view of Lemma 12731 and the discrete Poincare inequality 

implied by Lemma I5TT1 below we conclude that 

| < U,V >T | < C2|w|x|f U- 

In order to prove the positivity, we write < u,u >t in the form ( T2~8l 

using the algebraic inequality — 2ab ^ —a 2 — b 2 and the discrete boundary condition we obtain 

K&M a&£ K K&M a€£ K 

By the assumption (^3) one has that Vk,g ^ and we finally conclude that 

a&£ K 
< U,U >T^ 0. 



Next we recall a technical lemma presented in [T7], Lemma 8.2, which will be useful for the a 
priori estimates of the next section 

Lemma 2.5. Let B(s), s G R be defined by 

B(s) = P(s)s- [ 0{r)d 
Jo 



T. 



with f3 satisfying hypothesis (Hi). Then B(s) ^ T^ 2 ^' 
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3 A priori estimates. 



We define below an approximate solution of Problem (l)-(3). 
Definition 3.1. (Approximate solution) 

Let the sequence of{u n } G X^ , n G {1, . . . , N}, be a solution of the discrete problem ([74])- 
[17\). with 5t = T/N > 0. We say that the piecewise constant function u-r>,st : x [0, T] — ► R 
is an approximate solution of Problem (V) if 

Mx>,<5t(x, 0) = u° K for all x G K, 

u Vj6t (-K, t) = u\ for all (x, t) G K X (t„_i, £„]; 

i/i/e a/so define its approximate gradient by 

Vv,8t u vM x i *) = V©M n (x) for a// (x, t) e K x (t n _ u t n }. 

Lemma 3.1. f/\ priori estimate) Letuv,st be an approximate solution of Problem (l)-(3), 
then it is such that 

Ci = ^K)||i 2( n) +m(0)rM| o mm M F( M )| + ^||g||i 2(QT) ; 
moreover there exists C 2 > 0, stvc/j t/?at 

||/3(M©,5i)IU oo (0,T;L2(f7)) ^ C 2 . (30) 

Proof. Let m G [1, iV] be an arbitrary integer. Summing on n G {l,...,m} the equation 
([I5D with w n = w n for each n we obtain 



in 



E m w E <(£«) - /?(«r x )) + E u "> m " >^ + < m "> m " 

Jf£.M n=l n=l 

m rn 

+J2 6t E ™(K)u n K F{u n K ) = e E 



n=l XeA4 n=l ifeX 

Next, we consider the function 5 from Lemma I2~5l defined by 



B(u) = p(u)u - / p(r)dT. 
One can see that the following relation holds 



o 



B(u n K ) - BK- 1 ) = u n Mu n K ) - - / ifi(r) - ftu^dr 



n-l 
K 
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and since is nondecreasing we have that 

r X (/3(r)-/3«- 1 ))dr^O, 
JuT 1 

which implies 

m 

E m ^( B K) - = E m W E( 5 K) " 

ife.M -f£TeA4 n=l 

m 

< E m(K)E^(/3K)-/3K^))- 

KeM n=l 

In view of Lemma IZ51 we have that 

\§u* < B(„) < «0(„) < 

which yields 



^ll«D,ft(-.*m)lli»(n)-|ll)8(«o)llz>(n) < E m(K)J2^W(n n K )-/3(u n £ 1 )). 

^ If i- A A m — 1 



KeM n=l 

We remark that in view of the hypothesis (H5) one has 

uF(u) > M min F(u), 

since min < 0. The last statement of Lemma EH implies < u n ,u n >t^ 0. By the 

equation ( 123"|) and ("H 2 ) we finally conclude that 

^PWuvMi'^mjWhin) + MVv,stUv,8t\\h(nx(o,km)) 

m (31) 

^ + E E tt™(K)u n K q%, 

n=l if e.M 

where 

C = ^H/SMII^n) +m(fi)TM| o mm f 
Applying Cauchy-Schwarz and Young's inequality to the last term in (13"!]) leads to 

g^ll^ftC^OIIiafn) + All Vi?, jt'Uu.tft 1 1 1^ (nx (0,fcm)) ^ C+ ||wi>,<5t|U 2 (QT)ll9lU 2 (QT) 
< C + |ll«27,«lli»(0 T ) + 7^IMl! 2 (Q T )- 



We then obtain 

^||M©,5t|lL°°(0,T;L 2 (n)) ^ C + ^ ^ II U V,5t II L« (0,T;L 2 (fi)) + ^ II 9 II L 2 (Q T ) 



^ ~" '|2 ^ /-y i £ rr>|i . 1 1 2 ' n n'2 
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and ^ 

M\^T>M u T>,5t\\h(Q T ) ^ C + 2 T \\ UT> ' St \\L oo (0,T;L !l (Q)) + T^T II <?ll L 2 (Q T )- 

We now choose e = (3/{2T), which gives 

-^P\\ u V,&t\\ 2 L^{0,T;L^(Q,)) ^ & a nd ^ II ^V,StUv,5t || L 2 (Q T ) ^ C> 

where 

G = ^\\P(uo)\\h m + m{n)TM\ o min M F( M )| + ^||?||i w 

In order to prove the estimate on ||/3(Mx>,5t)lk°°(o,T;L 2 (n)) we sp I it (3 into a bounded and a 
Lipschitz continuous part by setting /3 = ft + ft, where 



and 



«'H« ^ AM-U, ^ (32) 



/3(P)-P(-P) p(P) + 0(-P) 

V{*) = { 2P S+ 2 |SKP (33) 

\s\ > P. 



We finally define 



Pi = Pi-y and ft = ft + 2/ . (34) 



we then remark that /3 = ft + ft; we remark that ft and ft are continuous and that ft is 
bounded by 2Cp, while ft is Lipschitz continuous with Lipschitz constant = max (ft (ftP) — 
P(-P))/2P). Which implies the L°°(0, T; L 2 (fi)) estimate 

||/3(«lV*)IU° o (0,T;L 2 (n)) ^ Wl3l(uT>,8t)\\L°°(0,T;L 2 (n)) 
+ ||ft(«Z>,5i) - ft(0)||z/~(0,T;L 2 (C)) + || ft (0) || L°° (0,T;L 2 (ft)) , 

so that 

||/3(^,<5t)||L°°(0,T;L 2 (n)) < 2/72(^)5(7^ + L^|| Uj), ffl || L oo (0,T;L 2 (n)) + m(fi) ^ | ft (0) | . 

Remark 3.1. (Extended discrete problem) Let s > and w G H M (VL) (sf. Definition 
\2.5\) . we consider the following extended one step problem. Find u E Xx> t o such that for all 

v e X Vfi : 

E, r , N /3(«k) — P{wk) 
m{K)VK r h < t>,-u >i? +S < V,U >t 
ot 

K€M (35) 
+s 22 m {K)v K F(u K ) = s 22 m (K) v K(lK- 
KeM KeM 
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It can be shown that the solution of the extended problem {35]) is bounded in the norm \-\x- 
More precisely, it is such that 

\mv,stUv,6t\\h(Q T ) ^K^II/3(^)lli 2(fi) +m(0)TM| o mm M F( M )| + ^|| g ||i 2(QT) ). (36) 

Theorem 3.1. (Existence of a discrete solution) The problem fl4] ) -(fT7] ) has at least one 
solution. 

Proof. Let (e^i^c^x-D ) be a family elements of Xx>fl, which components are defined by 
(ei)j = 5ij, where Sij is the Kronecker symbol. The system of nonlinear equations (ED-([l2D 
may be written in the form 

E{P{u n ) - Piu 71 - 1 )) + Au n + Cu n + 5tEF{u n ) = Q n , (37) 

where 

(i) u n , u n ~ l G X Vfi ; 

(ii) E is the diagonal matrix of the size card(Ai) + card(E int ) with elements 

(E) KjK = m(K) and (E)^ = 

for all K G M, a G £ int ; 

(iii) (3 and F are continuous mappings from Xp to itself naturally defined by 

(P(u))i = 0(ui) and (F(«)) i =F(« i ); 

(iv) ^4 and C are the diffusion matrix and the convection matrix respectively, with components 

Aij — 6t < ej, ej >^ and = 5t < ej, >r, 

(v) Q n G o is the source term, given by 

Q n K = 5tm(K)q n K and Q™ = 

for all if G .M, a G £j„ t ; 

Due to the coercivity of the bilinear form corresponding to the diffusion the matrix A is 
invertible; hence (I3T1) is equivalent to 

u n + A-\E{f3{u n ) - (3{u n ' 1 )) + Cu n + 5tEF{u n ) - Q n ) = 0. 

As it has been done in ( I3"5j) . we introduce the extended formulation 

u n + - Piu"- 1 )) + Cu n + 5tEF{u n ) - Q n ) = 0, (38) 
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with s G [0,1]. Moreover for a given u n 1 we define a continuous mapping H n : [0,l]xIp iO -> 
X V fi by 

H n (s,u) = sA-\E(f3(u) - /3(M n_1 )) + Cu + 6tEF(u) - Q n ). 

Then the equation (1H5|) can be written in the form u n + H n (s,u n ) = 0. In view of Remark 
13.11 the estimate (13"0~|) and the Lemma 12721 we have that 

5t\u\ 2 x < C, 

with some positive constant C, which does not depend on s. Setting R = sfCfbt + 1 we 
deduce that 

|u|x < R for all (s, u) G [0, 1] x X v>0 such that it + H n (s, u) = 0. 

Therefore the equation u + H n (s, u) = has no solutions on the boundary of the ball of 
radius R for s G [0, 1]. Next, we denote by d(Id + H n (s, •), Br, 0) the topological degree of 
the application Id + H n (s, ■) with respect to the ball Br and right-hand side 0. In view of 
the homotopy invariance of the topological degree and thanks to the fact that H n (0,u) = 
for all u G Xx>fi we have that 

d(Id + H n (s, •), B R , 0) = d(Id + H n {0, •), B R , 0) = 1 for all s G [0, 1], 

where we have applied [[11], Theorem 3.1 (dl) and (d3)]. Thus, by [[11], Theorem 3.1 (d4)], 
there exists u n such that u n + H n (l, u n ) = 0, so that u n is a solution of (|3"7|) . 



Theorem 3.2. (Uniqueness of the discrete solution) The solution of the problem (fl4])- 
( 07]) is unique. 

Proof. We give a proof by contradiction. Let ux>,st and u h k be two different solutions of ( TT4|) - 
(T17]) . such that w m = u m for all m = 1, . . . , n - 1, but u n ^ u n . We define r n = u n - u n . In 
view of (TT5]) with t> = r n we have that 

y m(f)r «MhM! + < r « r « > F + < r » r " > T 



+ £ m(i^K(F«) - = 0. 



We apply Lemma l2~4l as well as the assumptions ("Hi) and ("H 5 ) in order to estimate each 
term in the above equation. We obtain that 

W/St-F) ]T m(K)(r n K ) 2 + a\r n \ 2 x ^0, 
KeM 

where a is the coercivity constant. Finally, in view of the assumption ([5]) on the time step we 
deduce that 

\r n \x = 0- 
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4 Estimate on time translates 



To begin with we give two technical lemmas which will be useful for proving the estimate on 
time translates 

Lemma 4.1. Let T > 0, r G (0, T), N G N \ {0}, St = T/N be given and (a™)neN\{o} 

be a 

family of non negative real values. Let \s] denotes the smallest integer larger or equal to s. 
Then 

, T -r N 
\t/St]+l^n^[(t+T)/Si] n=l 

Proof. One has that 

»T-t rT-T rT-T 



I 



j J2 andt ^ [ J2 aHdt = / J2 am+ldt 

^° \t/St]+l^n^\(t+T)/8t\ ^° t/5t+l^n<(t+T)/St+l ^° tsim5t<t+T 



We remark that if \t/8t] + 1 > \(t + r)/St\, then the above inequality seal holds, with the 
left hand side term equal to zero. We define a characteristic function x( n ,ti,h) by 



X{n,t 1 ,t 2 ) 



1 if h ^ n5t < t 2 , 
otherwise. 



Then we obtain that 

cT-T AT— 1 rT-T N 



I 



_ ^ ± pi —T 



\t/5t\+l^n^\(t+T)/8t\ m=l " u m=l 

Lemma 4.2. Let T > 0, r e (0,T), iV G N \ {0}, <ft = T/N, ( G [0,r] be given and 
(a n )«eN\{o} ^ e a family of nonnegative real values. Let \s] denotes the smallest integer larger 
or equal to s. Then 



I 



a^Wdt^rjha". 

" \t/Si]+l^n^[(t+T)/St] n=l 



Proof. As in the proof of the previous Lemma we have that 

cT-t N rT-r 



I 



■ -^vi± , ^ i Jo 



rt/< 5 *i+ is s™^r( t +' r )/ 5t i n=i 

A simple change of variable implies 

N pT _ T N 



E/ 7 a^ t+ ^x(n,t,t + r)dt <W a^xfa s _ C , s _ C + r ) ds 

n=l ^ n=l ^ 

J2 am J2 s - C, s - C + r)ds. 

m=l n=l J (rn-l)St 
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In order to conclude the proof we remark that xi n i £2) = x{ n + m,ti + m8t, £2 + m8t) for 
all n, m G Z, £i,£2 G ^- which in turn implies that 

N /-m<5t N t-~n5t 

V/ - C,s ~ ( + r)ds = E/ x(-m,s-(,s-( + T)ds 

n=1 J(m-X)8t n=1 J~(n+l)St 



Theorem 4.1. Let V be a discretization ofVL in the sense of Definition 12.21 and let {ux>,st} 
be a solution of the discrete problem in sense of Definition I3.il Let also 9 ^ 9x> be given. 
Then there exists a positive constant C only depending on 9 such that 

/ / (uv,st(x, t + r) — u V jt(x, i)) 2 dxdt ^ Cr, (39) 
Jo Jn 

for all re (0,T). 

Proof. To begin with we use the hypothesis ("Hi) to obtain 

P / (uv,st( x > t + r) - u VjSt (x, t)) 2 dxdt 

Jo Jn 

= p[ T ~ T E m(K)(v}t T)m -u^fdt 

< f T E ™{K){uT T)m -4 /<5tl )(^(4 t+r)/M1 )-/3(4 /5tl ))d£ 

^° KeM 

= [ T ' T E HK){uT T)m - E - / w 1 ))* 

For a given /c and for all real t and r we define the following set 

n(£, r) = {n G N, [£/6£] + 1 ^ n < \(t + r) /5t] }, 
which can be empty. Then, the discrete equation ([7]) implies 

§_ r r i (u VM (x,t+T)- UVM (x,t)) 2 d X dt < r r e (4* +T)/ati 

Jo Jn Jo Xe-M 

■ H^{K)ql-Yi F KA^) + yK,^)-m{K)F{u n K ))dt. 

n£n(t,r) ct&£k 

Let us define the expressions Aoc, A fi and by 
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A D ,c = 



f T E st E ^ t+T)m - uT') E (*w« B ) + 

^° nen(t,r) KeM ae£ K 



rT 

A R = 



A S =f' T e ^E-w(^ +r)/5<i --r 1 )^, 

^° n£n(t,r) #eA4 

which we will estimate below. In view of ( TI3j) . (TIB]) and ( I17j) we obtain 

A D , C = /"^ ]T (jt(< u r(*+r)/«i _ u r*/«i , > F + < u r(*-Hr)/«i - u r*/«i , M « >t) ^ 

^° nen(t,r) 

In view of Lemma 12.41 we have that \ < u,v > F \ + \ < u,v > F | ^ C|w|x|^|x f° r a " u , v £ 
X V fl and since 2a6 $C a 2 + b 2 one has 

\A DtC \^C T T M(\u mT)m \x + \u {tm \x)\u n \xdt 

^° n6n(t,r) 

^° n€n(t,T) ^° nen(t,T) ^° nen(t,r) 

It follows from the estimate (l2~9"|) and the Lemmas 14711 and l4~2l that 

iV 
n=l 

Next, we estimate the term Ar, we remark that for all u,v G M it holds 

uF(u) < Lj?|v||«| < ^v 2 + ^L 2 F u 2 if it < 0, 

vF{u) ^ |f | max \F(u)\ ^-v 2 + - max F 2 (u) if < u ^ M, 
vF(u) ^ \v\(\F(u) - F(M)\ + \F(M)\) ^ \v\(L F \u -M\ + F(M)) 
^v 2 + ^L 2 F \u\ 2 + ^(L F M + F(M)) 2 \fu>M. 



Hence, 



E ™( k >kF(u k ) ^ \\v\\ 2 LHn) + h 2 F \\u\\ 2 LHn) + l -C F 
KeM 
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for all v, u e H M , where C F = -m(Q)( max F 2 (u) + (L F M + F(M)f). We obtain 

/T — T 
St (hv,st(; r*/^l)lli>(n) + ll«2>,«0, r(* + ^)/^l)llxa(n))d* 
nen(t,T) 

^ 5t(L||| M ^(-,t n )lli 2(n) + C F )dt. 

neti((,r) 

One more time it follows from the estimate (j25} and the Lemmas I47H and l4~2~l that 

1^1 ^r^5t(C|| U2?i5t (-,t n )||i 2(n) + C F ) ^rC 

n=l 

In the same way we proceed for the term \As\, one has that 

fT — T ft n f 

A S= EE// ( U K +T)m - uT')^ S) dxdsdt 



n£n(t,T) KeM 

and 

pT—t r t 



\a s \ < / x: E /" / ic^^+c^^n + c^*)) 2 ^* 

= / E |**(ll«23,*t(-, r*/^l)llz^(n) + r(*H-T)/**l)lli=»(n))^ 

^° nen(t,r) 
pT—T pt n p 

+ / EE/ / (<?( x > s )) 2 rfxM - 



n€n(t,T) KeM Jtn ~ 1 JA " 

In view of Lemmas l4~2l and l4~T1 we obtain 

Finally we use an a priori estimate (|29|) and the hypothesis (He) to conclude the proof. 

5 Estimate on space translates 

In this section we prove an estimate on the L 2 — norm of differences of space translates of the 
discrete solution. We state without proof two results from [14J, which are useful in our study. 



Lemma 5.1. Let d ^ 1, 1 ^ p < oo and f2 be an open bounded connected subset ofM. d . Let 
V be a mesh ofVt in the sense of Definition ^. 21 Let rj > be such that r\ $C dK,a/di,u ^ V 7 ? 
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for all a G M. a = {K, L}. Then, there exists q > p only depending on p and there exist a 
positive constant C , only depending on d, Vt, p and i] such that: 



for all u G H^i(Q). We recall that H_m(Q) C L 2 (f2) is the set of piecewise constant functions 
on the control volumes of the mesh. 

Lemma 5.2. Let d ^ 1 and Vt be a polyhedral open bounded connected subset ofW 1 . Let 
V = (Ai,£,V) be a discretization of £1 in the sense of Definition 12.21 and let u G Hm{&>)- 
Then, with notation of Definition ^. 51 ' 



where u is defined on the whole M. d , taking u = outside £1. 
Next we show that a similar inequality holds in every L p -norm. 

Lemma 5.3. Let d^l, l^p<co and f2 be an open bounded connected subset ofW 1 
and T > 0. Let V be a discretization of in the sense of Definition \2.2[ Let 77 > such 
that r] ^ dK,a/di,u ^ l/ 7 / f° r a ^ e = {K, L}. There exist C > and p > 0, only 
depending on d, p, £1 and 77 such that 

IK" + y) - w||lp(r«») < C\y\ p \\u\\ hpM} 
where u is defined on R. d , taking u = outside £1. 

Proof. In view of Lemma I5~H there exist q > p and a positive constant C such that 



(40) 



u(- + y) -u\\ L i {Rd) ^ |y|Vd|| 



u \\l,l,Mi 



(41) 



We apply the Interpolation Inequality [[lj, Theorem 2.11, p. 27] 

IK- + y) - «l|LP(Rrf) < IK- + y) - «ll£i (R <i)IK- + y) - u 



(42) 



i-p 

Li(R d )' 



(43) 




IK + y) - ^11 LP(R d ) < c\y\ p (\\u\\ 1)1M ) p (\\u\\ ltPM ) 

Applying Holder inequality we obtain that 

IK|l,l,./W ^ C II ^|| l,p,A4 

for some positive constant C . Then 

\\u{- + y) - u\\ LP{Rd) ^ C|y| p ||M||i iPiA 4. 
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Theorem 5.1. Let {T>h} be a family of discretizations in the sense of Definition \2.2\ and 
let 9 be a positive constant such that 9x> ^ 9 for all V E XV Let {u-r>,st} be a family of 
approximate solutions corresponding to V h and 5t = T/N for some N E N \ {0}. Then 
{uv,st} /s relatively compact in L 2 (Q T ). 

Proof. To begin with, we extend u v ^t by zero outside of Q T . Applying the Lemma [5751 with 
p = 2 yields 

\\uv,st(- +y,t) - ^x>,st(-^)llL 2 (R d ) < C\y\ p \\u VM (-,t)\\ 1:2M 
for some positive constants p > and C > 0. Integrating on (0,T) we obtain 

N 

\\uv,5t{- + y, ■) - UV,6t\\L2(R<ix(0,T)) ^ C \y\ 2p 5t\\uy^t(-,t n ) \\i,2,M- 

n=l 

Then in view of the Iemmas l2~3"l 12.21 and the estimate ( T2"9~j) we obtain the bound 

||Wp,Jt(- + y) - UT),8t\\L 2 (M. d x(0,T)) ^ C|y| P 7 

which, combined with (139]) gives 

||wr>,<5t(" + y ; ' + r ) — u V,8t\\L 2 {R d x{0,T)) 
< ||«23,5i(- +y, ■ +T) - Wx>,<Jt(- + y, -)l|L 2 (K d x(Q,T)) + II^D,<5t(- + Y, ~ ^Z>,<5< || L 2 (IR d x (0,T)) 

^C(vT + |y|'). 

Then the Frechet-Kolmogorov Compactness Theorem implies that the family {ux>,st\ is rela- 
tively compact in L 2 (R d x (0,T)) and thus in L 2 (Q T ). 

6 Convergence 

Theorem 6.1. Let {V h } be a family of discretizations in the sense of Definition \2.2\ and 
let 9 be a positive constant such that 9x> ^ 6 for all V E Vh. Let {ut>,8t} be a family of 
approximate solutions corresponding to Vh and St = T/N for some N E N \ {0}. Then 
there exists a subsequence of{ux> : 8t}, which we denote again by {ux>,st], sucn that ux>,st u 
strongly in L 2 (Q T ) as hv,5t — > 0, where u is a weak solution of Problem (V). Moreover 
u E L 2 (0, T; iif- (SI)) and Wv,8t Ur DM weakly converge in L 2 (Q T ) d to Vu. In the case that 
F is nondecreasing, the whole sequence {ux>,st} converges to the unique weak solution u of 
Problem (P). 

Proof. By Theorem 15.11 there exist a subsequence of {uv,st} that we still denote by {ux>,st} 
and a function u E L 2 (Q T ) such that Ur>,st u strongly in L 2 (Q T ) as hx>,5t — >■ (and 
also in L 2 (R d x (0,T)) taking u VM = outside of Q x (0,T)). In view of Q55]) there exists 
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a function G £ L 2 (Qx) d such that Vv,8tUv,8t weakly converge in L 2 (Qx) d to G along a 
subsequence as hx>,5t — > 0. In order to show that G = Vw we consider an arbitrary vector 
function w £ C([0, T]; C c °°(M d )) and the term T G defined by 

t g= / Vx> )( 5tMc,5t(x, *) • w(x, t) dxdt. 

JO JR d 

Using the Definition O and ([20]) we obtain that T<i = T<| + T<|, with 



N 

t g = E E E - <) n *,. ■ ™k 

n=l KeM <?eE K 



and 



t g = E E E RkA^)^ ■ I I w(x,t) dxdt, 

n=l KeM <reE K Jtn-i JD K ^ a 

where w^- = — - / / w(x, t) dxdt. We compare T G with T G defined by 

otm(K) J tn _ x J K 

N 

t g = E E E 5t ^)« - • w» 

1 /■*» /■ 

where w™ = — — / / w(x, t) d'jdt. One can see that 

- < E E E 5J P^« - <) 2 E E E *m*)^k& - <\ 



n=l Ke.M o-6£k ^ K ' a n=l KeM aeS K 

which leads to T G — )■ Tq as h v — > 0. Then, 

t g = ~ E E E 8tm{a)u n K n K , a -w n c = - [ [ u VM (x, t) V ■ w(x, t) dxdt. 

n=l KeM ae£ K 



JR d 



and we conclude that 



lim T G = — / / u(x, t)V ■ w(x, i) rfxcit. 



Next, we show that lim T G = 0. Thanks to (}2"2"j) we have that 

t g = E E E v(«>x,^ ■ / / ( w ( x - *) - w ^) rfx ^ 



n=l KeM ae£i< 
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Since w is a regular function, there exist a positive constant C = C(w) such that 

(w(x,t) - w£) dxdt\ <: C5t m ^l dK ' a {h v + 5t). 
t n _i « 

On the other hand by (1211) and in view of regularity of the mesh, we have that 



n n ti tl 

{R K . a u n f < 2d((^-^) 2 + |V^"| 2 | X(J r X ^ | 2 ) < 2rf((^^) 2 + 2 |V^| 2 )- 

«if,cr «K> dK,a 

Applying the Cauchy-Schwarz inequality inequality, we get 

lim Tp = 

This implies that the function G G L 2 (IR d x (0,T)) d is a.e. equal to Vu in M ' x (0,T). Since 
u = outside of Q, it follows that w G L 2 (0, T; #q(^))- 

Next we show that u is a weak solution of the problem (V). For this purpose, we introduce 
the function space 

$ = {y?e C 2,1 (n x [0,T]), ^ = 0on5Ox [0,T], y>(-,T) = 0}. 

Taking an arbitrary ip E we define the sequence of elements of Xt>,o 

cp n = P v cp{-, t n ) for alln E {1, . . . , N} 

which implies ip\ = y?(x^,t n ) and y?" = ^(x ff ,t n ). Next setting 

v n = for all n e {1,...,N} 

in (TI5]) . we obtain, also in view of ( TIB]) and (TIT]) that 



where 



t t = E E - ^ uT k 1 Wk 1 

n=l ATGA4 
A 7 " 

^ = E 5t E E - 

n=l KG.M o-g£k 
TV 

n=i /ce» o-efx 

TV 

^ = E 5t E ™(K) V n K l F{i 



and 



TV 



^ = E 5t E m{K)^ql. 

n=l K&M 

We successively search for the limit of each of these terms as h-p and k tend to zero. 



24 



6.1 Time evolution term 



Using discrete integration by parts and the fact that </?(x,T) = we obtain 



A' 



n=l ATeA4 ifeA4 

First we show that 

Y,rn{K)P{u K )ip° K ^ I /3( Mo (x))^(x,0) dx. 
KeM ^ n 

For this purpose we define 

T°=J2 m{K)P(u° K )<p° K - I /3(n (x))^(x,0) dx. 

Next we subtract J n (3(u° K )(p(x, 0) from each term to deduce that, 

T t= E / /3«)(^-^x,0))rfx- £ f((3(u (x))-(3(u° K )Mx,0)dx. (44) 

In view of the regularity of the test function ip £ C 2,l (Vt x [0, T\) we have that 

|^-^(x,0)| ^ Ch v for all x £ K 

and 

|p(x,0)| <C. 

We also remark that by ([6]) one has that \u° K \ ^ ||wo|U°°(n) and moreover the monotonicity 
hypothesis (Hi) implies that \(3(u^)\ ^ /3(||wo|U°°(n)) f° r a" K £ .M; consequently the first 
term on the right-hand side of ( 144|) tends to zero as hx> — > and the second term can be 
estimated by 

C E / I^K(x)) dx = C f |/3(u (x)) -/3(« 27l «(x,0))| dx. 

By the discrete initial condition ([5]) and the Definition 13. H one has 

/ |w (x) - uv,st{x-, 0) | dx ->• as /in ->■ 0, 

or in other words wd^O) converges strongly to w in L 1 (fi) as h-p — >■ 0. Hence a subsequence 
of {«B,ji(x, 0)}, which we still denote by {^ ^(x, 0)} converges to u (x) for a.e. x £ and 
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also (3(u v ,8t(x, 0)) -)> /3(« (x)) for a.e. x G fi. Since /3(uD 5 5t( x , 0)) < ||/3(tio(x))|| io0 (n) the 
Lebesgue dominated convergence theorem implies 

/ |/3(m (x)) - f3(u VM (x, 0))| dx^Oas hj> -> 0. 
Jq 

Thus T® — > as h-p — > 0. Next we prove that 

TV t 

E E - -»> / / /3(u(x,t))^(x,t) dxdt (45) 



n=l ATeJV/f 

as Zip and 5t — >■ 0. We define 



= E E m (*0(^r - V n K l )P{u n K ) - I ! /3(«(x,t))^(x,t) dxdt, 

n=l if G.M ^° ^ n 

and we add and subtract / / /3(w^)^(x, t) dxdt in each term to obtain 

Jtn-! JK 

T t = E E m (K)f3(u n K ) / (^(x^t) -<^(x,t)) dxdt 

n=l KeM Jt n -\ (45) 



+ / / (l3(u V; 6t(x, t)) - 0(it(x, t)))^ t (x, t) dxdt. 

Jo </!T2 

We have that for all x G X and all if G .M it holds 

|y^(x^,t) - <p*(x,t)| < C(h v ) 

where C(hx>) — >■ as /i£> — > 0. The absolute value of the first term on the right-hand side of 
(146]) is bounded by 

TV TV 
n=l K£M n=l KcM 

^C(^)Tm(fi) 1 / 2 ||/3(n^)|| LO o (0iT;L2(n)) , 

which tends to zero as h-p — > in view of the a priori estimate (13Uj) . Further, since |y?t(x, t)| ^ 
C^, we can estimate the absolute value of the second term in ( 146]) by 

C v f I |/3(«u f i t (x,t)) -/3(«(x,t))| dxdt <C^/ / |^i(u23,5t(x,t)) -/3i(u(x,t))| dxdt 

J0 JH in 

+0, / / |Aj(iiu,ft(x,t)) -/3 2 («(x,t))| dxdt, 
Jo Jq 

where /3i and (3 2 are given by (!32|) - (l34|) . Since Wd^j — >■ u strongly in L 2 (Qt), a subsequence 
of {uv,st}, which we still denote by {uv,st} converges to u a.e. in Q. The first term on 
right-hand side of the expression above converges to zero by the Lebesgue dominated con- 
vergence theorem. The convergence to zero of the second term can be deduced from the 
Lipschitz continuity of /3 2 an d the strong convergence of Ur>,st to u in L 2 (Q T ). 
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6.2 Convection term 



Next, we show that 

JV 



E = X> E E (^" 1 -^r 1 )^,^, (T + / / M (x,t)V(x)-V^(x,t) rfxrft -> (47) 

n=l KeMa££ K ^° ^ Q 



as ht>, St tend to zero. As it was done in in the proof of the Lemma l2~4l we write the left-hand 
side part of ( 147|) as 



TV JV 

E^EE fa* 1 - vt'Wk^ = E 5t E E v^r 1 - 

n=l ifS.M o-ef/f n=l i^GM cr££ K /, Q \ 

-E 5t E E v^r l -?r 1 )(«ir-«2), 

n=l JfeA4 o-e£K,Vft:, CT <0 

and the following estimate holds 

JV JV 

I E 5t E E tV^F 1 - vT^Wk - <)\ < k> • l|v|U~ C n> E^i^ui^u 

n=l JfeX o-efjf ,Vr:, ct <0 n=i 

JV 

^ Ch v ■ \\V\\ Looin) J2St\\Vvf n ' 1 \\LHn ) \\V D u n \\ L 2 {n) 

n=l 

JV 

< C7tx> ■ ||V||//x>(fi) 2j ^||Vl><^ n_1 ||i2(Q) + C/Vp ■ ||V||L°o(n)||VD,5fW2J )( Sf||i2(Q T ). 

71=1 

The second term in the right-hand side of the expression above is bounded be because of the 
a priori estimate (129|) and the first term can be controlled via the consistency of the discrete 
gradient given by Lemma EH] and the regularity of ip; indeed 

N ftn f N ftn f 

E / / l V ^ n_1 ( x )| 2 d * dt < 3 E / / IV^^x) - V^(x,t n _ 1 )| 2 d*dt 

n=l ''tn—l n=l ^ 

+3^/" / |Vy(x,t n _i) - V^(x,t)| 2 dxdt + 3 / / |V<^(x,t)| 2 rfxdt 

^ (C^t) + C^)Tm(fi) + 3|| V^||i 2(QT) ^ C, 

where C^St) tends to zero as St — > 0. Thus the second term in the right hand side of ( 145|) 
tends to zero as hp, St — > 0. Let us define Ei and £7 2 

JV jv . 

E i = E 5t E E tWtfsr 1 - + E 5t E «* / v^vo • v( X ) 

n=l ATG.M crg^A' n=l ATeM 
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£ 2 = E E w * I V^(x, t n -i) ■ V(x) dx.- I f u(x, £)V</?(x, t) ■ V(x) rfxdt. 

n=l KeM Jk J ° Jn 

so that also in view of (ffBj) lim^ r x t ->.n E = lim hT>i st->-oEi — \im hVj st^oE 2 . We will successively 
establish that E 1 and E 2 converges to zero as hx>,St — > 0. To begin with let us remark that 
integrating by parts in the expression of E\ yields E\ = E u — E 12 , where 



JV N 



n=l KeMae£ K n=l KeM " K 



re£ K 

and where 

N N 



E n = X> E E VkoVT^k - E> E E M * / </>0Mn-i)V(x) • n^dr. 



We first prove that lim E n = 0. 

h-p,8t-K} 



N 

En = E 6t E < / (^sr 1 - ^ x > v • v w) rfx - 

n=l iteA4 J K 



KeM 

in view of regularity of the function ip we obtain 

N 



|£n| ^C^E 5 * E K| / |V-V(x)| dx^Cyi / |^(x,t)V- V(x)| rfx 
n=l Ke.M ^ K J Qt 

Finally applying the Cauchy-Schwarz inequality yields 

l^nK ^||^|| l2(Qt) || v • v|| i2(QT) 

so that | .En | — > as hx> — > 0. Next we consider the term Ei 2 , which can be written as 

- f 
Ei2 = J2 5t J2 E u k / - ^(x,t„-i))V(x) • n K>(T d r 

71=1 K€Mae£ K ^ u 

In order to show that £i 2 — >■ as hx>,5t — >• we first remark that since n^ )CT = — for 
any pair of neighbor volumes K, L, and in view of the boundary condition on ip one has that 

E 5t E E ^rt=o 

7i=l ifeX o-e£x 



and also 

N 



X> E E < / ^(x, vi)v(x) • = o. 



n=i AreA'f ae£ K 
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Hence, the term E19. can be written as 



'12 
N 



^ — 1 A A ^ <7 



n=l KeM aeS K 

Therefore, in view of the regularity of cp and V we have that 

A' 



\E 12 \ <: CmaxV(x) ■ h^5t m(a 

n=l KeMae£ K 



Applying Cauchy-Schwarz inequality we obtain 

ib 12 k cdWvw.Mi: E E «m W <^^<£)i.(E E E 

n=l AG.A/! o-g£k if,<T n=l K&M a<=£ K ' 

In view of Lemma I2~2l we obtain 

\E 12 \ ^ CdmaxV(x)m(0)^T5 • V^w^lU^), 

so that in view of the a priori estimate (129]) one has that |£ , i 2 | — > as /i© — > 0, so that 
-Ei = -En — -E12 as hv, St — >■ 0. It remains to prove that E 2 converges to zero. 

Adding and subtracting / / w^Vy^x, t) ■ V(x) dxdt from each term of E 2 , yields 

Jt„-i Jk 

E 2 = J2 E / / <(V^(x,t n _i) - Vp(x,t)) • V(x) d*dt 

^ /■*" /" 

"EE/ / M x > *) - u ^^(x, t))Vp(x, £) ■ V(x) dxcft. 



n=l A'ex 



Finally, in view of the regularity of the a priori estimate (l29l . and to the fact that uv,st — > u 
strongly in L 2 (Qt) we conclude that \E 2 \ tends to zero as hx>,St — > 0. 



6.3 Diffusion term 

We show below that 

t d = E 5t E E (vsr 1 - ^r 1 )^^) - / / v ^(x, • A( X ) v M (x, t) d*dt 

n=l KeMae£ K ^° ^ Q 

tends to zero as hx>, St — > 0. In view of ( 12~3|) one has that 

- /■*» /" 

= E / / (V^ n_1 • A(x) V v u n - V^(x, t) • A(x)Vm(x, t)) efacd*. 



71=1 
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Adding and subtracting the term / V<^(x, t) ■ A(x)V©u n dx. we set Tp in the form 

Jn 

1% + Tp with 

- r tn r 

T d = H / ( V »^ n_1 - V ^( X > *)) • A(x) V»M n dxdt 
n= i J t n -i J n 

and 



Td = J2 f n f V^(x, t) • A(x)(V©M n - Vm(x, t)) dxdt. 

The term T|, tends to zero as hv,5t — > 0, since Vv,&t u v,8t tends to Vm weakly in L 2 (Q T ). 
On the other hand the term T|> can be written in the form T|, = Tp + T|, with 

n= l tn— 1 J H 

and 



T 5 D = J2 [ I V (^( x > *n-0 - <^( x > *)) ■ A(x) V©u n rfxdt. 
n — i J t n —i j n 

It follows from ( T2"§j) , Lemma 12.11 and the regularity of cp that and T|, tends to zero as 
h v , St -> and so do and T£. 



6.4 Reaction term 

Let us show that 



Tji — >■ [ j F(u(x,t))<p(x,t) dxdt 
Jo Jn 

as hv and k tend to zero. For this purpose, we introduce 

N r tn r 

t r = E E / / fofc" 1 - ^ ^ 



n=l if G.M 

iV 



T ^ = EE /" / ^x,*)(F«)-F(«(x,t)))dxdt. 

We obtain the convergence result similarly as for the time evolution term; more precisely we 
split the reaction F into a bounded and a Lipschitz continuous parts by setting 

F ( 8 ) = I F ^ < s ^ M F (s) = { ° ° ^ S ^ M 
1 otherwise, 2 1 F(s) otherwise, 
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and 

( \ ! < s «C M 

{ otherwise. 

We, then, define F\ = F\ — y and F 2 = F 2 + y which are both continuous; moreover is 
bounded by = max \F(s)\+F(M), and F 2 is Lipschitz continuous with Lipschitz constant 

Lp = max(Lp, F(M)/M). In view of the regularity of cp one has 

\(p(x,t) - v n K l \ < C ( h v + St) for all x G K,t G {t n -i,t n ], 

so that 

|T^| ^ C{h v + St) I [ iF^uvA^t)) + F 2 (u VM (*,t))\ dxdt 
Jo Jn 

< C{h v + 5t){C p Tm{Q) + L p T^m{n) l 2\\u VM \\mQ T )), 

which by the a priori estimate ( 121?]) implies that \T^\ — > as hp, St — > 0. Since </? is bounded 
we can estimate the second term as 

1^11 <C [ [ \F{uv,st{x,t)) - F{u{x,t))\ dxdt 
Jo Jn 

^ [ [ C\F 1 {u- D>S t{x,t))-F 1 {u{x,t))\dxdt+ [ f C\F 2 {u v , 5t {x,t))-F 2 {u{x,t))\ dxdt. 
Jo Jn Jo Jn 

The convergence is can be proved by applying similar arguments as for the time evolution 

term. 



6.5 Source term 

We deduce from the regularity of cp that 

T s = ^2^2 / <P(*K,tn~i)q{x,t) dxdt -> / / ip(x,t)q(x,t) dxdt 

as hv, St — > 0. 



6.6 Convergence to a weak solution of Problem (V) 

In view of Theorem I6.1l {n-n a/} strongly converges to u in L 2 (Q T ), with u G L 2 (0, T; Hq(CI)), 
and it follows from (!30|) that 0(u) G L°°(0, T; L 2 (fi)). Moreover we deduce from the density 
of the set $ in the set {<p G L 2 (0, T; /f^(fi)), cp t G L°°(Q T ),<p(-,T) = 0} that u is a weak 
solution of the continuous problem (V) in the sense of Definition 12.11 In the case that F is 
nondecreasing so that the solution of Problem (V) is unique (cf. Remark |2~T|) we conclude 
that the whole family {ux>,st} converges to u. 
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7 Numerical simulations 



In this section we present the results of numerical simulations. The purpose is to test our 
scheme in the case of problems with a known analytical solution. 



7.1 Numerical Test I 

We consider the equation 

d( M + M2 ) _ v • (A(x)Vu) + V • (V(x)u) + l -u^ = 

in the 3-dimensional space domain f2 = (0,2) x (0,1) x (0,1). We define the discontinuous 
A and V fields as follows 
For all %\ ^ 1 we set 

/I 0\ 
A = 1 and V = (4,0,0); 

\o o lj 

for all xi > 1 we set 




and V = (4,7,7). 

The initial and the Dirichlet boundary conditions are given by the exact solution 

u(x,t) = e^+^+^^. 

We remark that the velocity field V and the total flux A(x)Vu + V(x)« have a continuous 
normal trace across the discontinuity x = 1. We perform the simulations on 3-dimensional 
hexahedral meshes with random refinement (see Figured]), so that the mesh is nonmatching. 
In Table Q] we present simulation results with various mesh sizes h v and time steps k; we 
denote by Err the maximum relative error in L 2 -norm, namely 

„ IK,i(-,t„) - W(',*n) 11^(0) 

hrr = max — — . 

ne{l,...,N} \\u{-,t n )\\ L 2 {n) 

7.2 Numerical test II 

We consider a degenerate parabolic equation which possesses a traveling wave solution, namely 

- V • (5Vu) + V • {{v, 0, 0)u) = 
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Figure 1: Approximate solution on the nonmatching hexahedral mesh at t = 1 



N 


h 


# of elements 


# of faces 


Err 


50 


0.75 


165 


672 


0.03575 


100 


0.375 


837 


3324 


0.01432 


200 


0.1875 


3203 


11550 


0.00648 


400 


0.0938 


18533 


60633 


0.00305 



Table 1: Number of time steps N, mesh diameter hv, number of elements, number of faces 
and the relative error for nonmatching hexahedral meshes 

in the domain 

Q = (0, l) 3 and T = 1. 
This equation admits the following 1-dimensional exact solution 

u(x, y, t) = (1 - e ^ x - vt - p) f for x^vt+p, 

u(x, y, t) = for x > vt + p 

where p,v,8 are parameters still to be defined. We set p = 0.2, v = 0.8, and consider two 
values of 5, namely 0.01, 0.0001. The initial state is given by the exact solution at the time 
t — and we prescribe corresponding Dirichlet boundary conditions on the sides x = and 
x — 1. The null flux boundary condition is imposed on the remaining part of the boundary. 

Since the scheme does not preserve the maximum principle, it is necessary to define the function 
P(u) = for negative values as well, which leads us to set 0(—u) = —f3(u). Further one 
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has to solve the system of nonlinear equations 

m(K)((3(u n K ) - ^ul- 1 )) + k E F ^ n ) 

<re£ K 

+k Vk,o^ = St m{K)q n K , for all K G M, 

<tESk 

(^M + V K , a u n K(T ) + {F L>a {u n ) + Vl,„W2 = 0, for all a G E int 



(49) 



0, for all a £ £, 



ext- 



Since /3'(0) = +oo the Newton method can not be directly applied. In order to overcome this 
difficulty we introduce new discrete unknowns 

w n = /3(u n ), and thus u n = <p(w n ), where <p = 

In view of ( TlOl and (T25|) the nonlinear system becomes 

m{K){w n K - wl" 1 ) + St E A °°' 'ifPiO ~ 



E + = « ^ all i^GM 



inti 



E - + (^z>«) + ^(O) = 0, for all o" G £ 

<p(tu£) = 0, for all tr G £ ext . 

(50) 

The system (IHUjl depends on (w™) ae £ only through the terms ((p(w™)) ae £. This lead us to 
choose the discrete unknowns 

w n K = p(u") for all K G M and u n a for all cr G £, 

so that the system (1501 takes the form 

<)-*&)) 



<t,<t'&£k 



creSK 



E - <o + ( v kXk + y^o 



(51) 



+ E ^VK) - <') + (^VK) + ^X) = 0, for all a G 4 



inti 



u n a = 0, for all a E E. 



ext- 
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Figure 2: The approximate solution profiles at the time t = 0.5 for 5 = 0.01 and 5 = 0.0001 

Remark 7.1. We remark that the nonlinear system ( f4^) (or a linear one arising during the 
Newton 's procedure) has a special structure; more specifically, for each K G {1, . . . , card(Ai)} 
the K-th equation does not contain any unknown different from uk snd (u CT ) CTe g K (here we 
denote by K both the control volume and the index of the unknown uk ); therefore one can 
algebraically eliminate uk, so that the number of equations to solve becomes card(£). 

Since we do not impose many constraints on the mesh (in particular it can be nonconforming), 
it is not difficult to perform a local grid refinement. Finally note that there is a possibility 
to reduce the number of unknowns by using a method introduced in [14|; one can elimi- 
nate the interior interface unknowns (u a ) ae £ int by expressing them as a consistent barycentric 
combinations of the values uk- 
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